################################################################
#Replication code for the paper "Dominant party rule, 
#elections, and cabinet instability in African autocracies"
################################################################

library(dplyr)
library(broom)
library(ggplot2)
library(stargazer)
library(MASS)
library(multiwayvcov)
library(sandwich)
library(lmtest)
library(car)

#Load data directly from the web or enter file path for 
#Dataverse download

df_models <- read.csv(url('http://alexkroeger.weebly.com/uploads/5/6/1/7/56177689/bjpols_r2.csv'))

#df_models <- read.csv(your.file.path.here)

###############################################################
#Model formulas
###############################################################

fm1 <- as.formula(sacks ~ factor(gwf_party) + 
                    factor(gwf_military) +
                    factor(election_any_corrected) +
                    factor(opposition) + factor(attempt_corrected) +
                    sppop_L1  + factor(high_rents) + ln_cgdppc_L1 +
                    growth_L1 +
                    leader_tenure + I(leader_tenure^2) +
                    I(leader_tenure^3) + offset(log(cabinet_size_L1)))


fm1.5 <-  as.formula(sacks ~ factor(gwf_party)*factor(election_any_corrected) +
                       factor(gwf_military)*factor(election_any_corrected) +
                       factor(opposition) + factor(attempt_corrected) +
                       sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                       leader_tenure + I(leader_tenure^2) + 
                       I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fm2 <- as.formula(shuffled ~  factor(gwf_party) + 
                    factor(gwf_military) + 
                    factor(election_any_corrected) +
                    factor(opposition) + factor(attempt_corrected) +
                    sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                    growth_L1 +
                    leader_tenure + I(leader_tenure^2) + 
                    I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fm2.5 <- as.formula(shuffled ~ factor(gwf_party)*factor(election_any_corrected) +
                      factor(gwf_military)*factor(election_any_corrected) +
                      factor(opposition) + factor(attempt_corrected) +
                      sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                      growth_L1 +
                      leader_tenure + I(leader_tenure^2) + 
                      I(leader_tenure^3) + offset(log(cabinet_size_L1)))



#List of all formulas for lapply
allformulas <- list(fm1, fm1.5, fm2, fm2.5)

#######################################################
#Loop all model formulas for Table 1 through glm.nb
#######################################################

t1models <- lapply(allformulas, function(x){
  model <- glm.nb(x, data = df_models)
})

###########################################################
#Calculate clustered standard errors using multiwayvcov for
#table 1 models
###########################################################

t1SEs <- lapply(t1models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error
})

###########################################################
#Linear hypothesis tests for Model 1
###########################################################

linearHypothesis(t1models[[1]],
                 'factor(gwf_party)1 = factor(gwf_military)1',
                 vcov = cluster.vcov(t1models[[1]], df_models$obsid),
                 test = "F")


###########################################################
#Model comparisons and linear hypothesis tests for Model 2
###########################################################

#m1 vs m2
anova(t1models[[1]], t1models[[2]])


linearHypothesis(t1models[[2]], "factor(gwf_party)1 = factor(gwf_military)1",
                 vcov = cluster.vcov(t1models[[2]], df_models$obsid),
                 test = "F")



linearHypothesis(t1models[[2]],
                 "factor(gwf_party)1 + factor(election_any_corrected)1 + factor(gwf_party)1:factor(election_any_corrected)1 = factor(gwf_military)1 + factor(election_any_corrected)1 + factor(election_any_corrected)1:factor(gwf_military)1",
                 vcov = cluster.vcov(t1models[[2]], df_models$obsid),
                 test = 'F')


############################################################
#Linear hypothesis tests for Model 4
############################################################

linearHypothesis(t1models[[4]], "factor(gwf_party)1 = factor(gwf_military)1",
                 test = 'F', vcov = cluster.vcov(t1models[[4]], df_models$obsid))

linearHypothesis(t1models[[4]], "factor(gwf_party)1:factor(election_any_corrected)1 = factor(election_any_corrected)1:factor(gwf_military)1",
                 test = 'F', vcov = cluster.vcov(t1models[[4]], df_models$obsid))

########################################
#Table 1
########################################

table1 <- stargazer(t1models[1], t1models[2], t1models[3], t1models[4],
                    se = c(t1SEs[1], t1SEs[2], t1SEs[3], t1SEs[4]),
                    dep.var.labels = c("Dismissals","Horizontal Moves"),
                    order = c(1, 2, 3, 13, 14, 4, 7, 8, 9, 5, 6),
                    covariate.labels = c("Party Regime",
                                         "Military Regime",
                                         "Election",
                                         "Party X Election",
                                         "Military X Election",
                                         "Multiparty",
                                         "Resource Rents",
                                         "Ln(GDP per Capita)",
                                         "GDP Growth",
                                         "Coup Attempt",
                                         "Senior Partners",
                                         "Constant"),
                    omit = c("leader_tenure", "I(leader_tenure^2)", "I(leader_tenure^3)"),
                    column.sep.width = '.5pt',
                    no.space = T,
                    add.lines = list(c("Leader Tenure Polynomial?", "Yes", "Yes", "Yes", "Yes")))

#############################################################
#Figure 1 (Predictions from Models 2 and 4)
#############################################################


df_plot <- data.frame(gwf_party = c(1, 1, 0, 0, 0, 0),
                      gwf_military = c(0, 0, 0, 0, 1, 1),
                      election_any_corrected = c(0, 1, 1, 0, 0, 1),
                      opposition = c(1, 1, 1, 1, 1, 1),
                      attempt_corrected = c(0, 0, 0, 0, 0, 0),
                      sppop_L1 = c(mean(df_models$sppop_L1[df_models$gwf_party==1],),
                                   mean(df_models$sppop_L1[df_models$gwf_party==1],),
                                   mean(df_models$sppop_L1[df_models$gwf_personal==1],),
                                   mean(df_models$sppop_L1[df_models$gwf_personal==1],),
                                   mean(df_models$sppop_L1[df_models$gwf_military==1],),
                                   mean(df_models$sppop_L1[df_models$gwf_military==1],)),
                      high_rents = c(0,0,0,0,0,0),
                      ln_cgdppc_L1 = c(mean(df_models$ln_cgdppc_L1[df_models$gwf_party==1],),
                                       mean(df_models$ln_cgdppc_L1[df_models$gwf_party==1],),
                                       mean(df_models$ln_cgdppc_L1[df_models$gwf_personal==1],),
                                       mean(df_models$ln_cgdppc_L1[df_models$gwf_personal==1],),
                                       mean(df_models$ln_cgdppc_L1[df_models$gwf_military==1],),
                                       mean(df_models$ln_cgdppc_L1[df_models$gwf_military==1],)),
                      growth_L1 = c(mean(df_models$growth_L1[df_models$gwf_party==1],),
                                    mean(df_models$growth_L1[df_models$gwf_party==1],),
                                    mean(df_models$growth_L1[df_models$gwf_personal==1],),
                                    mean(df_models$growth_L1[df_models$gwf_personal==1],),
                                    mean(df_models$growth_L1[df_models$gwf_military==1],),
                                    mean(df_models$growth_L1[df_models$gwf_military==1],)),
                      leader_tenure = c(mean(df_models$leader_tenure[df_models$gwf_party==1],),
                                        mean(df_models$leader_tenure[df_models$gwf_party==1],),
                                        mean(df_models$leader_tenure[df_models$gwf_personal==1],),
                                        mean(df_models$leader_tenure[df_models$gwf_personal==1],),
                                        mean(df_models$leader_tenure[df_models$gwf_military==1],),
                                        mean(df_models$leader_tenure[df_models$gwf_military==1],)),
                      cabinet_size_L1 = c(mean(df_models$cabinet_size_L1[df_models$gwf_party==1],),
                                          mean(df_models$cabinet_size_L1[df_models$gwf_party==1],),
                                          mean(df_models$cabinet_size_L1[df_models$gwf_personal==1],),
                                          mean(df_models$cabinet_size_L1[df_models$gwf_personal==1],),
                                          mean(df_models$cabinet_size_L1[df_models$gwf_military==1],),
                                          mean(df_models$cabinet_size_L1[df_models$gwf_military==1],)))


pred1 <- data.frame(predict(t1models[[2]], newdata = df_plot, type = "response", se.fit = T, interval = 'confidence'))


pred2 <- data.frame(predict(t1models[[4]], newdata = df_plot, type = "response", se.fit = T))


pred1 <- pred1 %>%
  mutate(party = c('Party','Party','Personal','Personal', 'Military', 'Military'),
         election = c('No','Yes','Yes','No', 'No', 'Yes'))
pred2 <- pred2 %>%
  mutate(party = c('Party','Party','Personal','Personal', 'Military', 'Military'),
         election = c('No','Yes','Yes','No', 'No', 'Yes'))

pred1 <- pred1 %>%
  mutate(ul = fit + 1.96*se.fit, 
         ll = fit - 1.96*se.fit)

pred2 <- pred2 %>%
  mutate(ul = fit + 1.96*se.fit, 
         ll = fit - 1.96*se.fit)


plot1 <- ggplot(data = pred1, aes(x=as.factor(election), y=fit,
                                  ymin=ll, ymax=ul,group=party)) + 
  geom_linerange() + geom_pointrange() + 
  geom_hline(yintercept = 0, linetype = 'dashed', colour = 'black') + 
  facet_wrap(~party) +
  theme_bw() + xlab("Election Year?") + ylab("Dismissals") +
  scale_y_continuous(limits=c(0,11), breaks=c(0,1,2,3,4,5,6,7,8,9,10,11)) +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.minor = element_blank())

plot1



plot2 <- ggplot(data = pred2, aes(x=as.factor(election), y=fit,
                                  ymin=ll, ymax=ul,group=party)) + 
  geom_linerange() + geom_pointrange() + 
  geom_hline(yintercept = 0, linetype = 'dashed', colour = 'black') + 
  facet_wrap(~party) +
  theme_bw() + xlab("Election Year?") + ylab("Horizontal Moves") +
  scale_y_continuous(limits=c(0,11), breaks=c(0,1,2,3,4,5,6,7,8,9,10,11)) +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.minor = element_blank())

plot2



#############################################################################
#Plot average marginal effects for each covariate (Calculated in Stata)
#############################################################################

#Dismissals

marg1 <- read.delim(url('http://alexkroeger.weebly.com/uploads/5/6/1/7/56177689/marg2.txt'),
                    sep = ',', header = T)

#marg1 <- read.delim('your.file.path.here',
#                   sep = ",", header =T)

covorder <- c(1:10)

marg1 <- cbind(marg1, covorder)

marg1 <- arrange(marg1, covorder)

marg1 <- marg1 %>%
  mutate(ul = mfx+(1.96*se),
         ll = mfx-(1.96*se))

pmarg1 <- ggplot(marg1, aes(x = covorder, y = mfx,  ymin = ll,
                            ymax = ul)) +
  geom_pointrange() + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(limits = c(.5, 9.5),
                     breaks = c(1:10), 
                     labels = c("Party Regime", "Military Regime", "Election", "Multiparty",
                                "Coup Attempt","Senior Partners", "High Rents",
                                "Ln(GDP per Capita)", "GDP Growth", 
                                "Leader Tenure")) +
  xlab("") +
  ylab("") +
  scale_y_continuous(limits = c(-4.1, 4.5),
                     breaks = c(-4, -3, -2, -1, 0, 1, 2, 3, 4)) +
  coord_flip() +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.minor = element_blank())

pmarg1



#Horizontal Reshuffles

marg2 <- read.delim(url('http://alexkroeger.weebly.com/uploads/5/6/1/7/56177689/marg4.txt'),
                    sep = ',', header = T)



#marg2 <- read.delim('your.file.path.here',
#                    sep = ",")

marg2 <- cbind(marg2, covorder)

marg2 <- arrange(marg2, covorder)

marg2 <- marg2 %>%
  mutate(ul = mfx + (1.96*se),
         ll = mfx -(1.96*se))

pmarg2 <- ggplot(marg2, aes(x = covorder, y = mfx,  ymin = ll,
                            ymax = ul)) +
  geom_pointrange() + 
  geom_hline(yintercept = 0, linetype = "dashed") +
  scale_x_continuous(limits = c(.5, 9.5),
                     breaks = c(1:10), 
                     labels = c("Party Regime", "Military Regime","Election", "Multiparty",
                                "Coup Attempt","Senior Partners", "High Rents",
                                "Ln(GDP per Capita)", "GDP Growth", 
                                "Leader Tenure")) +
  xlab("") +
  ylab("") +
  scale_y_continuous(limits = c(-4.1, 4.5),
                     breaks = c(-4, -3, -2, -1, 0, 1, 2, 3, 4)) +
  coord_flip() +
  theme_bw() +
  theme(axis.line = element_line(colour = "black"),
        panel.grid.minor = element_blank())

pmarg2

#############################################################
#Split Sample Fixed Effects Models
#############################################################


df_party <- df_models %>%
  filter(gwf_party==1)

df_per <- df_models %>%
  filter(gwf_personal==1)

df_mil <- df_models %>%
  filter(gwf_military==1)

fmfe1 <- as.formula(sacks ~  election_any_corrected +
                      attempt_corrected + sppop_L1 +  
                      ln_cgdppc_L1 + growth_L1 +
                      leader_tenure + I(leader_tenure^2) +
                      I(leader_tenure^3) + as.factor(obsid) + 
                      offset(log(cabinet_size_L1)))


fmfe2 <- as.formula(shuffled ~  election_any_corrected +
                      attempt_corrected +  sppop_L1 +
                      ln_cgdppc_L1 + growth_L1 +
                      leader_tenure + I(leader_tenure^2) +
                      I(leader_tenure^3) + as.factor(obsid) +
                      offset(log(cabinet_size_L1)))

feformulas <- list(fmfe1, fmfe2)

################################################
#Party Regime Fixed Effects Models
################################################

t2partymodels <- lapply(feformulas, function(x){
  models <- glm.nb(x, data = df_party)
})

##################################################
#Party Regime Fixed Effects Model Standard Errors
##################################################

t2partySEs <- lapply(t2partymodels, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_party$obsid))
  cl <- tidy(cl)
  cl$std.error
})

###################################################
#Personalist Fixed Effects Models
###################################################

t2permodels <-lapply(feformulas, function(x){
  models <- glm.nb(x, data = df_per)
})

#########################################################
#Personalist Fixed Effects Model Standard Errors
#########################################################

t2perSEs <- lapply(t2permodels, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_per$obsid))
  cl <- tidy(cl)
  cl$std.error
})


###################################################
#Military Fixed Effects Models
###################################################

t2milmodels <-lapply(feformulas, function(x){
  models <- glm.nb(x, data = df_mil)
})

#########################################################
#Military Fixed Effects Model Standard Errors
#########################################################

t2milSEs <- lapply(t2milmodels, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_mil$obsid))
  cl <- tidy(cl)
  cl$std.error
})



##########################################################
#Table 2
##########################################################


table2 <- stargazer(t2partymodels[1],t2permodels[1],
                    t2milmodels[1], t2partymodels[2],
                    t2permodels[2], t2milmodels[2],
                    se = c(t2partySEs[1], t2perSEs[1],
                           t2milSEs[1], t2partySEs[2],
                           t2perSEs[2], t2milSEs[2]),
                    dep.var.labels = c("Dismissals","Horizontal Moves"),
                    order = c(1, 2, 3, 4, 5),
                    covariate.labels = c("Election",
                                         "Coup Attempt",
                                         "Senior Partners",
                                         "Ln(GDP per Capita)",
                                         "GDP Growth", 
                                         "Constant"),
                    omit = c("leader_tenure", "obsid"),
                    column.sep.width = '2pt',
                    add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))



###################################################################
#Random Effects Negative Binomial Regression
###################################################################
library(lme4)


#Center variables to aid convergence

dfme <- df_models %>%
  mutate(size_cen = cabinet_size_L1 - mean(cabinet_size_L1),
         tenyrs_cen = leader_tenure - mean(leader_tenure),
         gdppc_cen = ln_cgdppc_L1 - mean(ln_cgdppc_L1),
         sppop_cen = sppop_L1 - mean(sppop_L1))


me1 <- glmer.nb(sacks ~ gwf_party*election_any_corrected +  
                  opposition  + attempt_corrected + sppop_L1 + 
                  high_rents + gdppc_cen + tenyrs_cen +
                  offset(log(cabinet_size_L1)) + (1|obsid), 
                data = dfme,
                control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 10000000)))
summary(me1)

me2 <- glmer.nb(sacks ~ gwf_party*election_any_corrected + 
                  gwf_military*election_any_corrected +
                  opposition  + attempt_corrected + sppop_L1 + 
                  high_rents + gdppc_cen + tenyrs_cen +
                  offset(log(cabinet_size_L1)) + (1|obsid), 
                data = dfme,
                control=glmerControl(optimizer="bobyqa", optCtrl = list(maxfun = 10000000)))
summary(me2)



me3 <- glmer.nb(shuffled ~ gwf_party*election_any_corrected +  
                  opposition + attempt_corrected + sppop_L1 +
                  high_rents + gdppc_cen + tenyrs_cen + 
                  offset(log(cabinet_size_L1)) + (1|obsid), 
                data = dfme,
                control=glmerControl(optimizer="nloptwrap", optCtrl = list(maxfun = 1000000000)))
summary(me3)


me4 <- glmer.nb(shuffled ~ gwf_party*election_any_corrected +
                  gwf_military*election_any_corrected +
                  opposition + attempt_corrected + sppop_L1 +
                  high_rents + gdppc_cen + tenyrs_cen + 
                  offset(log(cabinet_size_L1)) + (1|obsid), 
                data = dfme,
                control=glmerControl(optimizer="nloptwrap", optCtrl = list(maxfun = 10000000000)))

summary(me4)

#############################################################
#Random Effects Model Tables
#############################################################


tableA1 <- stargazer(me1, me2, me3, me4,
                     dep.var.labels = c("Dismissals", "Horizontal Moves"),
                     order = c(1, 2, 10, 3, 11, 4, 5, 6, 7, 8, 9),
                     covariate.labels = c("Party Regime", 
                                          "Election",
                                          "Party X Election",
                                          "Military",
                                          "Mil. X Election",
                                          "Multiparty",
                                          "Coup Attempt",
                                          "Senior Partners",
                                          "Resource Rents",
                                          "Ln(GDP per Capita)",
                                          "Leader Tenure", 
                                          "Constant"),
                     column.sep.width = '2pt',
                     no.space = TRUE)

#############################################################
#Recode Party Regimes and Control for Polity2
#############################################################

#Code pure party and party-hybrid regimes

hybrids <- df_models %>%
  mutate(pure_party = ifelse(gwf_regimetype=='party-based',1,0),
         hybrid_party = ifelse(gwf_party==1 & pure_party==0, 1, 0),
         party_personal = ifelse(gwf_regimetype=='party-personal', 1, 0),
         party_military = ifelse(gwf_regimetype=='party-military', 1, 0),
         military_personal = ifelse(gwf_regimetype=='military-personal', 1, 0),
         military = ifelse(gwf_military==1 & military_personal != 1, 1, 0))



fmhybrid1 <-  as.formula(sacks ~ factor(pure_party)*factor(election_any_corrected) +
                           factor(hybrid_party)*factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fmhybrid2 <-  as.formula(sacks ~ factor(pure_party)*factor(election_any_corrected) +
                           factor(hybrid_party)*factor(election_any_corrected) +
                           factor(gwf_military)*factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fmhybrid3 <- as.formula(shuffled ~ factor(pure_party)*factor(election_any_corrected) +
                          factor(hybrid_party)*factor(election_any_corrected) +
                          factor(opposition) + factor(attempt_corrected) +
                          sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                          growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))


fmhybrid4 <- as.formula(shuffled ~ factor(pure_party)*factor(election_any_corrected) +
                          factor(hybrid_party)*factor(election_any_corrected) +
                          factor(gwf_military)*factor(election_any_corrected) +
                          factor(opposition) + factor(attempt_corrected) +
                          sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                          growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fmhybrid5 <-  as.formula(sacks ~ factor(pure_party)*factor(election_any_corrected) +
                           factor(party_personal)*factor(election_any_corrected) +
                           factor(party_military)*factor(election_any_corrected) +
                           factor(military_personal)*factor(election_any_corrected) +
                           factor(military)*factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))


fmhybrid6 <-  as.formula(shuffled ~ factor(pure_party)*factor(election_any_corrected) +
                           factor(party_personal)*factor(election_any_corrected) +
                           factor(party_military)*factor(election_any_corrected) +
                           factor(military)*factor(election_any_corrected) +
                           factor(military_personal)*factor(election_any_corrected) +
                           factor(military)*factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))



hybridparty_formulas <- list(fmhybrid1, fmhybrid2, fmhybrid3, fmhybrid4, fmhybrid5, fmhybrid6)

#Loop formulas for pure party models
hybridparty_models <- lapply(hybridparty_formulas, function(x){
  models <- glm.nb(x, data = hybrids)
})

#Loop clustered standard errors for pure party models

hybridparty_SEs <- lapply(hybridparty_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, hybrids$obsid))
  cl <- tidy(cl)
  cl$std.error  
})


#Pure party table

pure_hybrid_party <- stargazer(hybridparty_models[1], hybridparty_models[2],
                               hybridparty_models[3], hybridparty_models[4],
                               se = c(hybridparty_SEs[1], hybridparty_SEs[2],
                                      hybridparty_SEs[3], hybridparty_SEs[4]),
                               order = c(1, 3, 4, 2, 14, 15, 16, 5, 6, 7, 8, 9, 10),
                               dep.var.labels = c("Dismissals", "Horizontal Moves"),
                               covariate.labels = c("Pure Party",
                                                    "Hybrid Party",
                                                    "Military",
                                                    "Election",
                                                    "Pure Party X Election",
                                                    "Hybrid Party X Election",
                                                    "Military X Election",
                                                    "Multiparty",
                                                    "Coup Attempt",
                                                    "Senior Partners",
                                                    "Resource Rents",
                                                    "Ln(GDP per Capita)",
                                                    "Growth",
                                                    "Constant"),
                               omit = "leader_tenure",
                               no.space = TRUE,
                               column.sep.width = '2pt',
                               add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))





#Control for Polity2

fmpolity1 <- as.formula(sacks ~ gwf_party*election_any_corrected +
                          gwf_military*election_any_corrected +
                          opposition + polity2 + attempt_corrected + 
                          sppop_L1  +
                          high_rents + ln_cgdppc_L1 + growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))




fmpolity2 <- as.formula(shuffled ~ gwf_party*election_any_corrected +
                          gwf_military*election_any_corrected +
                          opposition + polity2 + attempt_corrected + 
                          sppop_L1 +
                          high_rents + ln_cgdppc_L1 + growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))

#Polity formulas
polity_formulas <- list(fmpolity1, fmpolity2)

#Loop polity formulas through glm.nb
polity_models <- lapply(polity_formulas, function(x){
  models <- glm.nb(x, data = df_models)
})

#Loop clustered standard errors for polity models
polity_SEs <- lapply(polity_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error   
})



polity <- stargazer(polity_models[1], polity_models[2],
                    se = c(polity_SEs[1], polity_SEs[2]),
                    order = c(1, 3, 2, 14, 15, 4, 5, 6, 7, 8, 9, 10),
                    dep.var.labels = c("Dismissals", "Horizontal Moves"),
                    covariate.labels = c("Party Regime",
                                         "Military Regime",
                                         "Election",
                                         "Party X Election",
                                         "Military X Election",
                                         "Multiparty",
                                         "Polity2",
                                         "Coup Attempt",
                                         "Senior Partners",
                                         "Resource Rents",
                                         "Ln(GDP per Capita)",
                                         "Growth",
                                         "Constant"),
                    omit = "leader_tenure",
                    column.sep.width = '2pt',
                    no.space = TRUE,
                    add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))


#############################################################
#Table of Leaders for Supplementary Appendix
#############################################################


leaders <- df_models %>%
  group_by(obsid) %>%
  summarise(leader = first(leader),
            country = first(country),
            cowcode = first(cowcode),
            start = min(year), 
            end = max(year),
            regime = last(gwf_regimetype)) %>%
  arrange(country, start) %>%
  ungroup() %>%
  dplyr::select(-obsid)


countries <- leaders %>%
  group_by(cowcode) %>%
  summarise(country = first(country))

#Fix Capitalization of Country Names

maketitle = function(txt){
  theletters = strsplit(txt,'')[[1]]
  wh = c(1,which(theletters  == ' ') + 1)
  theletters[wh] = toupper(theletters[wh])
  paste(theletters,collapse='')
}

library(stringr)

titles <- str_trim(tolower(as.character(leaders$country)))

leaders$country <- sapply(titles, maketitle)

leaders$country <- gsub("[[:punct:]]", "", leaders$country)

leaders$country <- gsub("kinshasa", "-Kinshasa", leaders$country)

leaders$country <- gsub(" brazzaville", "-Brazzaville", leaders$country)

leaders$country <- gsub("bissau", "-Bissau", leaders$country)

library(xtable)

print(xtable(leaders), include.rownames = F)


#############################################################
#Summary Statistics for Supplementary Appendix
#############################################################

stargazer(df_models,
          covariate.labels = c("Dismissals",
                               "Dismissals Reduced",
                               "Horizontal Shuffles",
                               "Horizontal Reduced",
                               "Promotions",
                               "Demotions",
                               "Party Regime",
                               "Personal Regime",
                               "Military Regime",
                               "Election",
                               "Coup Attempt",
                               "Cabinet Size",
                               "Senior Partners",
                               "Ln(GDP per Capita)",
                               "Resource Rents",
                               "GDP Growth",
                               "Leader Tenure",
                               "Polity 2",
                               "Multiparty",
                               "Parliamentary"),
          omit = c("cowcode", "obsid", "year",
                   "I(leader_tenure^2)",
                   "leader_tenure2",
                   "leader_tenure3",
                   "I(leader_tenure^3)",
                   "N",
                   "gwf_duration",
                   "not_first",
                   "party_notfirst",
                   "pres",
                   'partytype',
                   'party_first'))



##############################################################
#Controlling for Presidential vs. Parliamentary Systems
##############################################################


pres_parl_sacks <- as.formula(sacks ~ gwf_party*election_any_corrected +
                                gwf_military*election_any_corrected +
                                opposition + parl + attempt_corrected + 
                                sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
                                leader_tenure + I(leader_tenure^2) + 
                                I(leader_tenure^3) + offset(log(cabinet_size_L1)))



pres_parl_shuffles <- as.formula(shuffled ~ gwf_party*election_any_corrected +
                                   gwf_military*election_any_corrected + 
                                   opposition + parl + attempt_corrected +  
                                   sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
                                   leader_tenure + I(leader_tenure^2) + 
                                   I(leader_tenure^3)+ offset(log(cabinet_size_L1)))


#Presidential/Parliamentary system formulas
presparl_formulas <- list(pres_parl_sacks,
                          pres_parl_shuffles)


#Loop Pres/Parl formulas through glm.nb

presparl_models <- lapply(presparl_formulas, function(x){
  models <- glm.nb(x, data = df_models)
})

#Loop clustered standard errors for pres/parl models

presparl_SEs <- lapply(presparl_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error
})

#Pres/Parl Table


pres_parl <- stargazer(presparl_models[1], presparl_models[2], 
                       se = c(presparl_SEs[1], presparl_SEs[2]),
                       order = c(1, 3, 2, 14, 15, 4, 5, 6, 7, 8, 9, 10),
                       dep.var.labels = c("Dismissals", "Horizontal Moves"),
                       covariate.labels = c("Party Regime",
                                            "Military Regime",
                                            "Election",
                                            "Party X Election",
                                            "Military X Election",
                                            "Multiparty",
                                            "Parliamentary",
                                            "Coup Attempt",
                                            "Senior Partners",
                                            "Resource Rents",
                                            "Ln(GDP per Capita)",
                                            "GDP Growth",
                                            "Constant"),
                       omit = "leader_tenure",
                       column.sep.width = '2pt',
                       no.space = TRUE,
                       add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))

###############################################################
#Models on Reduced Sample of Ministers
##############################################################


reduced1 <-  as.formula(sacks_reduced ~ gwf_party*election_any_corrected +
                          gwf_military*election_any_corrected +
                          opposition + attempt_corrected + 
                          sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))

reduced2 <- as.formula(shuffled_reduced ~ gwf_party*election_any_corrected +
                         gwf_military*election_any_corrected +
                         opposition + attempt_corrected +
                         sppop_L1 + high_rents + ln_cgdppc_L1 + growth_L1 +
                         leader_tenure + I(leader_tenure^2) + 
                         I(leader_tenure^3) + offset(log(cabinet_size_L1)))


#Reduced minister sample formulas
reduced_formulas <- list(reduced1, reduced2)

#Loop Reduced minister formulas through glm.nb
reduced_models <- lapply(reduced_formulas, function(x){
  models <- glm.nb(x, data = df_models)
})


#Loop clustered standard errors

reduced_SEs <- lapply(reduced_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error
})


#Reduced minister sample table

reduced_table <- stargazer(reduced_models[1], reduced_models[2],
                           se = c(reduced_SEs[1], reduced_SEs[2]),
                           dep.var.labels = c("Dismissals","Horizontal Moves"),
                           order = c(1, 3, 2, 13, 14, 4, 5, 6, 7, 8, 9),
                           covariate.labels = c("Party Regime",
                                                "Military Regime",
                                                "Election",
                                                "Party X Election",
                                                "Military X Election",
                                                "Multiparty",
                                                "Coup Attempt",
                                                "Senior Partners",
                                                "Resource Rents",
                                                "Ln(GDP per Capita)",
                                                "GDP Growth",
                                                "Constant"),
                           omit = c("leader_tenure", "I(leader_tenure^2)", "I(leader_tenure^3)"),
                           column.sep.width = '2pt',
                           no.space = TRUE,
                           add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))



##################################################
#Promotions and Demotions
##################################################



fmpromote <-  as.formula(promotion ~ gwf_party*election_any_corrected +
                           gwf_military*election_any_corrected +
                           opposition + attempt_corrected +
                           sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))

fmdemote <-  as.formula(demotion ~ gwf_party*election_any_corrected +
                          gwf_military*election_any_corrected +
                          opposition + attempt_corrected + 
                          sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))

#Promotion and demotion formulas
prodemo_formulas <- list(fmpromote, fmdemote)

#Loop promotion and demotion formulas through glm.nb
prodemo_models <- lapply(prodemo_formulas, function(x){
  model <- glm.nb(x, data = df_models)
})

#Loop clustered standard errors for promotions and demotions
prodemo_SEs <- lapply(prodemo_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error  
})


#Promotions and demotions table


promote_demote <- stargazer(prodemo_models[1], prodemo_models[2],
                            se = c(prodemo_SEs[1], prodemo_SEs[2]),
                            dep.var.labels = c("Promotions","Demotions"),
                            order = c(1, 3, 2, 13, 14, 4, 5, 6, 7, 8, 9),
                            covariate.labels = c("Party Regime",
                                                 "Military Regime",
                                                 "Election",
                                                 "Party X Election",
                                                 "Military X Election",
                                                 "Multiparty",
                                                 "Coup Attempt",
                                                 "Senior Partners",
                                                 "Resource Rents",
                                                 "Ln(GDP per Capita)",
                                                 "GDP Growth",
                                                 "Constant"),
                            omit = c("leader_tenure", "I(leader_tenure^2)", "I(leader_tenure^3)"),
                            column.sep.width = '2pt',
                            no.space = T,
                            add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))

##################################################
#Predicted Demotions Plot
##################################################


pred3 <- data.frame(predict(demote, newdata = df_plot, type = 'response', se.fit = T))

pred3 <- pred3 %>%
  mutate(party = c('Party','Party','Personal/Military','Personal/Military'),
         election = c('No','Yes','Yes','No'))

pred3 <- pred3 %>%
  mutate(ul = fit + 1.96*se.fit, 
         ll = fit - 1.96*se.fit)

################################################################
#Founding and Nonfounding Party Regime Leaders
################################################################




fm2lead1 <-  as.formula(sacks ~ party_notfirst*election_any_corrected +
                          party_first*election_any_corrected +
                          gwf_military*election_any_corrected +
                          factor(opposition) + factor(attempt_corrected) +
                          sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                          leader_tenure + I(leader_tenure^2) + 
                          I(leader_tenure^3) + offset(log(cabinet_size_L1)))



fm2lead2 <- as.formula(shuffled ~ party_notfirst*election_any_corrected +
                         party_first*election_any_corrected +
                         gwf_military*election_any_corrected +
                         factor(opposition) + factor(attempt_corrected) +
                         sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                         growth_L1 +
                         leader_tenure + I(leader_tenure^2) + 
                         I(leader_tenure^3) + offset(log(cabinet_size_L1)))


#Founding/nonfounding leader formulas
nfl_formulas <- list(fm2lead1, fm2lead2)

#Loop fl/nfl formulas through glm.nb
nfl_models <- lapply(nfl_formulas, function(x){
  models <- glm.nb(x, data = df_models)
})

#Loop clustered standard errors for fl/nfl models
nfl_SEs <- lapply(nfl_models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error
})

#FL/NFL Table



table_fl_nfl <- stargazer(nfl_models[1], nfl_models[2],
                          se = c(nfl_SEs[1], nfl_SEs[2]),
                          dep.var.labels = c("Dismissals","Horizontal Moves"),
                          order = c(1, 3, 4, 2, 14, 15, 16, 5,
                                    6, 7, 8, 9, 10),
                          covariate.labels = c("Party (NFL) ",
                                               "Party (FL)",
                                               "Military",
                                               "Election",
                                               "Party (NFL) X Election",
                                               "Party (FL) X Election",
                                               "Military X Election",
                                               "Multiparty",
                                               "Coup Attempt",
                                               "Senior Partners",
                                               "Resource Rents",
                                               "Ln(GDP per Capita)",
                                               "GDP Growth",
                                               "Constant"),
                          omit = c("leader_tenure", "I(leader_tenure^2)", "I(leader_tenure^3)"),
                          column.sep.width = '2pt',
                          no.space = T,
                          add.lines = list(c("Cubic Splines?", "Yes", "Yes", "Yes", "Yes")))

#########################################################################
#Replicate Table 1 with cubic polynomial of regime tenure
#########################################################################



regimeten1 <- as.formula(sacks ~ factor(gwf_party) + 
                           factor(gwf_military) +
                           factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1  + factor(high_rents) + ln_cgdppc_L1 +
                           growth_L1 +
                           leader_tenure + I(leader_tenure^2) +
                           I(leader_tenure^3) + 
                           gwf_duration + I(gwf_duration^2) + 
                           I(gwf_duration^3) +
                           offset(log(cabinet_size_L1)))


regimeten2 <-  as.formula(sacks ~ factor(gwf_party)*factor(election_any_corrected) +
                            factor(gwf_military)*factor(election_any_corrected) +
                            factor(opposition) + factor(attempt_corrected) +
                            sppop_L1  + factor(high_rents) + ln_cgdppc_L1 + growth_L1 +
                            leader_tenure + I(leader_tenure^2) + 
                            I(leader_tenure^3) + 
                            gwf_duration + I(gwf_duration^2) + 
                            I(gwf_duration^3) +
                            offset(log(cabinet_size_L1)))



regimeten3 <- as.formula(shuffled ~  factor(gwf_party) + 
                           factor(gwf_military) +
                           factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                           growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) +
                           gwf_duration + I(gwf_duration^2) + 
                           I(gwf_duration^3) +
                           offset(log(cabinet_size_L1)))


regimeten4 <- as.formula(shuffled ~ factor(gwf_party)*factor(election_any_corrected) +
                           factor(gwf_military)*factor(election_any_corrected) +
                           factor(opposition) + factor(attempt_corrected) +
                           sppop_L1 + factor(high_rents) + ln_cgdppc_L1 +
                           growth_L1 +
                           leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) +
                           gwf_duration + I(gwf_duration^2) + 
                           I(gwf_duration^3) +
                           offset(log(cabinet_size_L1)))



#List of all formulas for lapply
allformulas.1 <- list(regimeten1, regimeten2, regimeten3, 
                      regimeten4)

#######################################################
#Loop all model formulas for Table 1 through glm.nb
#######################################################

t1.1models <- lapply(allformulas.1, function(x){
  model <- glm.nb(x, data = df_models)
})

###########################################################
#Calculate clustered standard errors using multiwayvcov for
#table 1 models
###########################################################

t1.1SEs <- lapply(t1.1models, function(x){
  cl <- coeftest(x, cluster.vcov(x, df_models$obsid))
  cl <- tidy(cl)
  cl$std.error
})


########################################
#Table 1
########################################

table1.1 <- stargazer(t1.1models[1], t1.1models[2],
                      t1.1models[3], t1.1models[4],
                      se = c(t1.1SEs[1], t1.1SEs[2],
                             t1.1SEs[3], t1.1SEs[4]),
                      dep.var.labels = c("Dismissals","Horizontal Moves"),
                      order = c(1, 3, 2, 16, 17, 4, 5, 6, 7, 8, 9),
                      covariate.labels = c("Party Regime",
                                           "Military Regime",
                                           "Election",
                                           "Party X Election",
                                           "Military X Election",
                                           "Multiparty",
                                           "Coup Attempt",
                                           "Senior Partners",
                                           "Resource Rents",
                                           "Ln(GDP per Capita)",
                                           "GDP Growth",
                                           "Constant"),
                      omit = c("leader_tenure", "I(leader_tenure^2)", "I(leader_tenure^3)",
                               "gwf_duration", "I(gwf_duration^2)", "I(gwf_duration^3)"),
                      column.sep.width = '2pt',
                      no.space = T,
                      add.lines = list(c("Cubic Splines, Leader?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                                       c("Cubic Splines, Regime", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")))

#############################################################
#Drop One Leader/Country Analyses
#############################################################



dflead <- df_models %>%
  group_by(obsid) %>%
  summarise(leader = first(leader),
            sacks = mean(sacks),
            shuffled = mean(shuffled),
            gwf_party = median(gwf_party))

dfcountry <- df %>%
  group_by(cowcode) %>%
  summarise()

mean(dflead$sacks[dflead$gwf_party==1],)


mean(dflead$sacks[dflead$gwf_party==0],)


#drop individual leaders and countries and recompute results



fmdrop1.2 <-  as.formula(sacks ~ gwf_party*election_any_corrected +
                           gwf_military*election_any_corrected +
                           opposition + attempt_corrected +
                           sppop_L1  + high_rents + ln_cgdppc_L1 + 
                           growth_L1 + leader_tenure + I(leader_tenure^2) + 
                           I(leader_tenure^3) + offset(log(cabinet_size_L1)))



leaders <- as.list(dflead$obsid)

countries <- as.list(dfcountry$cowcode)

drop1 <- function(droplist){
  
  reduced <- df %>%
    filter(obsid != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
}


rse1 <- function(droplist){
  
  reduced <- df %>%
    filter(obsid != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
  cl1 <- coeftest(model, cluster.vcov(model, reduced$obsid))
  cl1 <- tidy(cl1)
  rse <- cl1$p.value[2]
}

rse1_int <- function(droplist){
  
  reduced <- df %>%
    filter(obsid != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
  cl1 <- coeftest(model, cluster.vcov(model, reduced$obsid))
  cl1 <- tidy(cl1)
  rse <- cl1$p.value[14]
}

drop2 <- function(droplist){
  
  reduced <- df %>%
    filter(cowcode != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
}


rse2 <- function(droplist){
  
  reduced <- df %>%
    filter(cowcode != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
  cl1 <- coeftest(model, cluster.vcov(model, reduced$obsid))
  cl1 <- tidy(cl1)
  rse <- cl1$p.value[2]
}

rse2_int <- function(droplist){
  
  reduced <- df %>%
    filter(cowcode != droplist)
  
  model <- glm.nb(fmdrop1.2, data = reduced)
  
  cl1 <- coeftest(model, cluster.vcov(model, reduced$obsid))
  cl1 <- tidy(cl1)
  rse <- cl1$p.value[14]
}


all1 <- lapply(leaders, drop1)

dropone1 <- data.frame(sapply(all1, function(x) summary(x)$coefficients[2, 1:4]))

df_dropone1 <- data.frame(t(dropone1))

names(df_dropone1) <- c("estimate", "se", "z", "p")

cluster1 <- lapply(leaders, rse1)

cluster1 <- data.frame(cluster1)

cluster1 <- data.frame(t(cluster1))

names(cluster1) <- "clsep"

df_dropone1 <- data.frame(cbind(df_dropone1, cluster1))


pdo1e <- ggplot(data = df_dropone1, aes(x = estimate)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime Estimate") +
  ylab("Density")

pdo1e



pdo1p <- ggplot(data = df_dropone1, aes(x = clsep)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime Estimate P-Value") +
  ylab("Density")

pdo1p


#Party Regime x Election interaction

dropone1.5 <- data.frame(sapply(all1, function(x) summary(x)$coefficients[13, 1:4]))


df_dropone1.5 <- data.frame(t(dropone1.5))

names(df_dropone1.5) <- c("estimate", "se", "z", "p")

cluster1.5 <- lapply(leaders, rse1_int)

cluster1.5 <- data.frame(cluster1.5)

cluster1.5 <- data.frame(t(cluster1.5))

names(cluster1.5) <- "clsep"

df_dropone1.5 <- data.frame(cbind(df_dropone1.5, cluster1.5))

pdo1eint <- ggplot(data = df_dropone1.5, aes(x = estimate)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime X Election Estimate") +
  ylab("Density")

pdo1eint


pdo1pint <- ggplot(data = df_dropone1.5, aes(x = clsep)) +
  geom_density() +
  theme_bw() +
  xlab("Party X Election Estimate P-Value") +
  ylab("Density")

pdo1pint



##################################################
#Dropping Individual Countries
##################################################

#Party regime coefficient and p-value

all2 <- lapply(countries, drop2)

dropone2 <- data.frame(sapply(all2, function(x) summary(x)$coefficients[2, 1:4]))

df_dropone2 <- data.frame(t(dropone2))

names(df_dropone2) <- c("estimate", "se", "z", "p")


cluster2 <- lapply(countries, rse2)

cluster2 <- data.frame(cluster2)

cluster2 <- data.frame(t(cluster2))

names(cluster2) <- "clsep"

df_dropone2 <- data.frame(cbind(df_dropone2, cluster2))


pdo2e <- ggplot(data = df_dropone2, aes(x = estimate)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime Estimate") +
  ylab("Density")

pdo2e



pdo2p <- ggplot(data = df_dropone2, aes(x = clsep)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime Estimate P-Value") +
  ylab("Density")

pdo2p



#Party regime x election interaction estimate and p-value


dropone2.5 <- data.frame(sapply(all2, function(x) summary(x)$coefficients[13, 1:4]))

df_dropone2.5 <- data.frame(t(dropone2.5))

names(df_dropone2.5) <- c("estimate", "se", "z", "p")

cluster2.5 <- lapply(countries, rse2_int)

cluster2.5 <- data.frame(cluster2.5)

cluster2.5 <- data.frame(t(cluster2.5))

names(cluster2.5) <- "clsep"

df_dropone2.5 <- data.frame(cbind(df_dropone2.5, cluster2.5))



pdo2eint <- ggplot(data = df_dropone2.5, aes(x = estimate)) +
  geom_density() +
  theme_bw() +
  xlab("Party Regime X Election Estimate") +
  ylab("Density")

pdo2eint



pdo2pint <- ggplot(data = df_dropone2.5, aes(x = clsep)) +
  geom_density() +
  theme_bw() +
  xlab("Party X Election Estimate P-Value") +
  ylab("Density")

pdo2pint




#########################################################
#Extreme bounds analyses
#########################################################

library(ExtremeBounds)


cluster_ses <- function(nbmodel){
  vcov_mat <- cluster.vcov(nbmodel, df_models$obsid)
  ses <- sqrt(diag(vcov_mat))
  return(ses)
}

df_models$reg_type <- ifelse(df_models$gwf_party==1, 2, NA)

df_models$reg_type <- ifelse(df_models$gwf_military==1, 3, df_models$reg_type)

df_models$reg_type <- ifelse(df_models$gwf_personal==1, 1, df_models$reg_type)

m2eba <- eba(formula = sacks ~ factor(reg_type)*election_any_corrected| opposition + attempt_corrected +
               sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
               leader_tenure + leader_tenure2 + 
               leader_tenure3 +  cabinet_size_L1,
             data = df_models, se.fun = cluster_ses, reg.fun = glm.nb)


hist(m2eba, variables = c('factor(reg_type)2', 'factor(reg_type)3',
                          'election_any_corrected',
                          'factor(reg_type)2:election_any_corrected',
                          'factor(reg_type)3:election_any_corrected'),
     main = c('Party Regime', 
              'Military Regime',
              'Election',
              'Party X Election', 
              'Military x Election'))


print(m2eba)

m4eba <- eba(formula = shuffled ~ factor(reg_type)*election_any_corrected| opposition + attempt_corrected +
               sppop_L1  + high_rents + ln_cgdppc_L1 + growth_L1 +
               leader_tenure + leader_tenure2 + 
               leader_tenure3 +  cabinet_size_L1,
             data = df_models, se.fun = cluster_ses, reg.fun = glm.nb)

hist(m4eba, variables = c('factor(reg_type)2', 'factor(reg_type)3',
                          'election_any_corrected',
                          'factor(reg_type)2:election_any_corrected',
                          'factor(reg_type)3:election_any_corrected'),
     main = c('Party Regime', 
              'Military Regime',
              'Election',
              'Party X Election', 
              'Military x Election'))


print(m4eba)